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Q,^ , We have simulated the temporal evolution of pressure due to capillary and viscous forces in 

' two-phase drainage in porous media. We analyze our result in light of macroscopic flow equations 

for two-phase flow. We also investigate the effect of the trapped clusters on the pressure evolution 
and on the effective permeability of the system. We find that the capillary forces play an important 
role during the displacements for both fast and slow injection rates and both when the invading fluid 
is more or less viscous than the defending fluid. The simulations are based on a network simulator 
modeling two-phase drainage displacements on a two-dimensional lattice of tubes. 
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, I. INTRODUCTION 

i 

r]^ ■ Fluid flow in porous media such as sand, soil and fractured rock is an important process in nature and has a huge 
number of practical applications in engineering. Most often mentioned is oil recovery and hydrology. Fluid flow in 
^ I porous media has also been of great interest in modern physics. In particular, the different structures of the interface 

■ ^ ' between the fluids in two-phase displacements have been extensively studied. Despite this attention there are still 

many open questions concerning fluid flow in porous media. 
■ In this paper we report on simulations of the temporal evolution of pressure during two-phase drainage in a model 
porous medium, consisting of a two-dimensional lattice of tubes. The network model has been developed to measure 
the time dependence of different physical properties and to study the dynamics of the fluid movements. Especially, 
' we focus on the dynamics of the temporal evolution of the pressure due to capillary and viscous forces and the time 
^ i dependence of the front between the two liquids. The discussion is restricted to drainage displacement, i.e. the process 
OO where a non- wetting fluid displaces a wetting fluid in a porous medium. 

, During the last two decades an interplay between experimental results and numerical simulations has improved 
the understanding of the displacement process. It has been shown that the different structures observed when 
changing the physical parameters of the fluids like viscosity contrast, wettability, interfacial tension and displacement 
rate [Q-^ divide into three flow regimes. These three major regimes are referred to as viscous flngering |^,^ stable 
\ displacement and capillary fingering . There exist statistical models such as DLA |^ , anti-DLA ||^ and invasion 
. percolation pcf that reproduce the basic domains in viscous fingering, stable displacement and capillary fingering 
^ ■ respectively. However, these models do not contain any physical time for the front evolution and they cannot describe 

■ ^ , the crossover between the different fiow regimes. 

^sj" Much effort has gone into making better models whose properties are closer to that of real porous media. This 
f~i has resulted in several network simulators, modeling fluid flow on a lattice of pores and throats [^|J|,|l^-jl^ . Most 
O-i of the network models have been used to obtain new information on the different flow regimes and to study the 
^ statistical properties of the displacement structures. Some others, have been used to calculate macroscopic properties 
like fluid saturations and relative permeabilities and compare them with corresponding experimental data. However, 
to our knowledge no network model has been capable of simulating the dynamics of the fluid pressures during the 
displacement process. This is the goal of the present work. 

We have simulated the temporal evolution of the pressure in all the three regimes of interest: viscous fingering, 
stable displacement and capillary fingering. The injection rate in the displacements has been systematically varied 
and we have analyzed the behavior of the pressure in the crossover between the three regimes. Moreover, we discuss 
what effect trapped clusters have on the evolution of the pressure in the system and we relate the data to macroscopic 
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flow equations. We find that capillary forces play an important role in two-phase flow at both high and low injection 
rates. 

The paper is or gan ized as follows. In Sec. || we present the network model, in Sec. [II we show the simulation 
results and in Sec. IV we present our conclusions and related the findings to macroscopic quantities. 



II. THE NETWORK MODEL 

The network model has been presented in [20|] and we will restrict ourself to a short sketch here. 



A. Geometry and boundary conditions 

The porous medium is represented by a square lattice tubes connected at the nodes. At each node four tubes meet 
and there is no volume assigned to the nodes: the tubes represent the volume of both pores and throats. The tubes 
are cylindrical with length d. Each tube is assigned an average radius r which is chosen at random in the interval 
[Aid, X2d\, where < Ai < A2 < 1. The randomness of the radii represents the disorder in the model. 

The liquids flow from the bottom to the top of the lattice and we implement periodic boundary conditions in the 
horizontal direction. The pressure difference between the bottom row and the top row defines the pressure across the 
lattice. Gravity effects are neglected, and consequently we consider horizontal flow in a two-dimensional network of 
tubes. 



B. Fluid flow through the network 



Initially, the system is filled with a defending fiuid with viscosity /ii. The invading fluid with viscosity fi2 is injected 
along the bottom row with a constant injection rate. We model drainage, i.e. the invading fiuid is non-wetting and the 
defending fluid is wetting. Furthermore, we assume that the fluids are incompressible and immiscible. Consequently, 
the volume flux is conserved everywhere in the lattice and a well deflned interface develops between the two phases. 

The capillary pressure pc due to the interface between the non- wetting and wetting fiuid inside a tube (a meniscus) 
is given by the Young-Laplace law 

27 

Pc — COS0 . (1) 

r 

Here r is the radius of the tube, 7 is the interfacial tension and 6 denotes the wetting angle between the non-wetting 
and wetting phases. 9 is in the interval (0,7r/2) for drainage displacements. 

We assume that the tubes are hour glass shaped where the effective radii of the tubes follow a smooth function. 
Thus, the capillary pressure becomes a function of the position of the meniscus in the tube and we assume that the 
Young-Laplace law takes the form 

2t 

= — [1 - cos(27rx)] . (2) 
r 

Here x is the dimensionless value of the meniscus' position in the tube (0 < i < 1), and 9 = (perfect wetting). 
From Eq. (|^) = at the ends of the tube while Pc approaches its maximum of 47/r in the middle of the tube. 

The volume flux qij through a tube from the ith to the jth node in the lattice (Fig. ^ is found from the Washburn 
equation for capillary flow | |2l[ | . As an approximation we treat the tubes as if they were cylindrical and obtain 

nrf: kij 1 

Ii] = -li^Prj-Pc)- (3) 

Me// d 

Here kij is the permeability of the tube given by kij = i^fj/^ where is the average radius of the tube. Apij = pj —pi 
is the pressure difference between the ith and jth node. A tube partially filled with both of the liquids is allowed to 
contain either one or two menisci leading to four different arrangements as shown in Fig. The effective viscosity 
of those tubes, denoted as Me// in Eq- (H), becomes a sum of the amount of each fluid multiplied by their respective 
viscosities. The total capillary pressure, Pc in Eq. (^, is the sum of the capillary pressures of the menisci that are 
inside the tube. The absolute value of the capillary pressure of each meniscus is given by Eq. (0), while its sign 
depends on whether the meniscus is pointing upwards like in Fig. ^(a) or downwards like in Fig. §(b). In the simple 
case where the tube only contains one meniscus (Fig. ^ fJ-eff = fJ-2Xij + fj,i{l — Xij) and pc = Pc- For a tube without 
menisci Pc = 0, and Eq. (||) reduces to that describing Hagen-Poiseuille flow with /Xe// = fii or ^2- 
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C. Determining the flow field 



There is no volume assigned to the nodes giving conservation of volume flux at each node 

3 

The summation on j runs over the nearest neighbor nodes to the ith node while i runs over all nodes that do not 
belong to the top or bottom rows, that is, the internal nodes. 

Eqs. (||) and constitute a set of linear equations which are to be solved for the nodal pressures pj with the 
constraint that the pressures at the nodes belonging to the upper and lower rows are kept fixed. The set of equations 
is solved by using the Conjugate Gradient method [Q. 

We want to study the dynamics of the pressure fluctuations at constant displacement rate. Therefore, we need to 
find the pressure across the lattice for a desired injection rate and then use that pressure to solve fluid flow through 
the network. For two-phase displacement the pressure across the lattice AP is related to the injection rate Q through 
the equation 

Q = AAP + B . (5) 

Here A and B are parameters depending on the geometry of the medium and the current configuration of the liquids. 
The first part of Eq. (^) is simply Darcy's law for one phase flow, while the last part B results from the capillary 
pressure between the two phases. As long as the menisci do not move B is constant. 

The pressure Apy across each tube can be related to the pressure across the lattice AP. All the equations calculating 
the fluid flow in the system, have the functional form f{x) = ax + b. As a consequence Ap^ becomes a function of 
AP, 

Apy=r„AP + ny. (6) 

Tij is a dimensionless quantity depending on the mobilities (k/^ieff) of the tubes and fly is a function of the capillary 
pressures of the menisci inside the tubes. If no menisci are present Ily is zero. Eq. ^) can easily be deduced for two 
cylindrical tubes with different radii connected in series. 

By inserting Eq. (^) into Eq. (^) we obtain after some algebra a relation between the local flow rate qij and the 
pressure AP across the network 

Qij = a^jAP + bij . (7) 

The parameter iiij is proportional to and the mobility {kij / fi^ff) of the tube ij. bij contains the capillary pressures 
of the menisci. 

The solution due to a constant injection rate can now be summarized into the following steps: 

1. We first find the nodal pressures for two different pressures AP and AP applied across the lattice. 

2. From the two solutions of the nodal pressures the corresponding injection rates Q and Q and the local flow 
rate g', and g-'- are calculated. 

3. A and B is calculated by solving the two equations obtained when inserting AP , Q , AP and Q into Eq. (||). 

4. The pressure AP across the lattice for the desired Q is then calculated by using Eq. (^. 

5. This AP is inserted in Eq. (Q) to get the local flow qij. Note that the parameters ciij and bij must first be found 
by solving the two equations obtained by insterting q'^j , AP q'/j and AP from step 2 into Eq. (|^) . 



D. Moving the menisci 



A time step At is chosen such that every meniscus is allowed to travel at most a maximum step length Axmax 
during that time step. In each time step we check whether or not a meniscus crosses a node. If this happens, the 
time step is redefined such that this meniscus stops at the end of the tube. 

A meniscus reaching the end of a tube is moved into the neighbor tubes according to some defined rules. These 
rules take care of the different fiuid arrangements that can appear around the node. Basically, the non-wetting fluid 
can either invade into or retreat from the neighbor tubes as shown in Fig. ||(a) and (b) respectively. In Fig. ||(a) the 
non- wetting fluid approaches the node from below (drainage). When the meniscus has reached the end of the tube 
(position 1), it is removed and three new menisci are created at position 6 in the neighbor tubes (position 2). The 
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distance 5 is about 1-5% of the tube length d. The small distance 5 avoids that the created menisci at position 2 
immediately disappear and move back to the initial position 1 in tubes where the flow direction is opposite to the 
direction of the invading fluid. The total time lapse is adjusted to compensate the instantaneous change in local 
volume of the fluids when the menisci move a distance S. The time lapse is adjusted such that the total volume of 
the fluids always is conserved. 

Fig. §(b) shows the opposite case when the non- wetting fluid retreats into a single tube (imbibition). As Fig. |^(b) 
shows the properties of imbibition should not be neglected as long as the menisci can travel in both directions. 
Our approximation in Fig. ^(b) cannot handle important properties found in imbibition such as fllm flow and snap 
off |^,|2^. However, in drainage which is what we are focusing on, arrangement (b) will appear rarely compared to (a). 
For that reasion, any futher description of imbibition than the one presented in Fig. ^(b), does not seem necessary. 

Summarized, the procedure for each time step At is: 

1. The nodal pressures pj are determined. 

2. The pj^s are related to the desired injection rate Q from Eqs. (||) and (|^). 

3. The local flow rate in each tube is computed by using Eq. (^). 

4. The local flow rates are used to calculate the time step At such that only one meniscus reaches the end of a tube 
or travel at most the step length Axmax during that time step. 

5. The menisci are updated according to At. The total time lapse is recorded before the whole procedure is repeated 
for the new fluid configuration. 



III. SIMULATIONS 



In two-phase fluid displacement there are mainly three types of forces: viscous forces in the invading fluid, viscous 
forces in the defending fluid and capillary forces due to the interface between them. This leads to two dimensionless 
numbers that characterize the flow: the capillary number Ca and the viscosity ratio M. 

The capillary number describes the competition between capillary and viscous forces. It is defined as 

a = (8) 

where Q (cm^/s) denotes the injection rate, /z (Poise) is the maximum viscosity of the two fluids , E (cm^) is the 
cross section of the inlet and 7 (dyn/cm) is the interfacial tension between the two phases. E is calculated by taking 
the product of the length of the inlet and the mean thickness of the lattice due to the average radius of the tubes. 

M defines the ratio of the viscosities of the two fluids and is given by the invading viscosity /i2 divided with the 
defending viscosity /ii: 

M = ^ . (9) 

Ml 

In the simulations the pressure across the lattice AP is given by Eq. da) as 



Since B is due to the capillary pressure of the menisci, we define —B/A as the global capillary pressure of the system. 
Peg ■ Peg includes the menisci surrounding the trapped cluster of defending fluid (cluster menisci) as well as the menisci 
belonging to the front between the invading and defending fluid (front menisci). 

In addition to Peg we calculate Pcf, the capillary pressure averaged along the front. Pc/ consists only of the capillary 
pressures due to the front menisci and we deflne it as 



1 ^ 

^^/ = ]^ElPcl- (11) 



Here the index a addresses the tubes in the lattice and in the summation a runs over all tubes containing a meniscus 
that belong to the front. N is the number of such tubes, is the capillary pressure of the front meniscus in tube a. 
See Fig. ^ for an example of a structure obtained after a simulation with low capillary number. The tubes containing 
a front meniscus or belonging to the trapped cluster of defending fluid are identified by running a Hoshen-Kopelman 
algorithm ||23] on the lattice. 
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The displacement structure in Fig. ^ shows that the front is a comphcated path connecting the left and the right 
boundaries of the lattice. The width of the front w is found by taking the standard deviation of the distances between 
all the front tubes and the average position of the front. 

For every simulation we have calculated AP and Peg as functions of time. For some of the simulations we have also 
computed the average capillary pressure P^f along the front and analyzed the behavior of A in Eq. (j^). As will be 
discussed below, AP, Peg and Pcf are strongly correlated and A seems to obey surprisingly simple relations. 

We have performed drainage simulations in each of the regimes of interest: viscous fingering, stable displacement and 
capillary fingering with M = 1.0x10^"^, 1.0 and 1.0x10^ respectively. The injection rate was systematically varied 
for each of the viscosity ratios. The simulations were performed with parameters as close as possible to experiments 
performed in |2^] . The length d of all tubes in the lattices were set equal to 1 mm and the radii r of the tubes were 
chosen randomly in the interval 0.05d < r < d. The interfacial tension was set to 7 = 30dyn/cm and the viscosities 
of the defending and the invading fluids varied between 0.01 P (~ water) and 10 P (~ glycerol). 



A. Viscous fingering, M < 1 

Our study of this regime consists of a series of simulations with constant viscosity ratio M — 1.0x10^'^. The 
injection rates and the capillary number used in the different simulations are listed in Table ^ To save computation 
time, most of the simulations were performed on a lattice of 25 x 35 nodes. However, one was performed on a lattice 
of 60 X 80 nodes and the resulting structure of this simulation is shown in Fig. |[ 

In viscous fingering the principal force is due to the viscous forces in the defending fluid. The pattern formation in 
Fig. H shows that the invading fluid creates typical fingers into the defending fluid. The pressure across the lattice, 
AP, shown in Fig. H, decreases as the less viscous fluid invades the system. The slope of the average decreasing 
pressure function is non-trivial and results from the fractal development of the fingers. 

The global capillary pressure Peg, also shown in Fig. fluctuates around a mean value of about 1 x 10'^ dyn/cm'^. 
The fluctuation is strongly correlated in both time and amplitude to the noise in AP. As will be discussed below, the 
variations correspond to the changes in the capillary pressure as a meniscus invades into or retreats from a tube. 

We have calculated Pcf for every simulation performed on the lattices of 25 x 35 nodes. The result for four of the 
simulations is shown in Fig. ^ together with the global capillary pressure Peg. The pressures in Fig. ^ are normalized 
by dividing them with the average threshold pressure of the tubes. The average threshold pressure is defined as 27/ (r) 
where (r) is the mean radius of the tubes. The mean threshold pressure is about 1.1 x lO'^dyn/cm^ in all simulations, 
since the radii of the tubes always are chosen randomly in the interval [0.05, 1.0] mm. Note that in Fig. 0, Peg has been 
subtracted by 1000 dyn/cm^ before it was normalized. This avoids the pressure functions to overlap at low capillary 
numbers when they are plotted in the figure. 

Fig. ^ shows that the fluctuations of Pc/ are correlated in time to the fluctuations of Peg- At lowest Ca = 3.5 x 10~^, 
Pcf and Peg are even indistinguishable. For all simulations except that at lowest capillary number the amplitude of 
the fluctuations in Pcf decreases as time increases. For high injection rates Pc/ is found to approach 1, which is the 
mean threshold pressure. 

In the simulation at Ca — 3.5 x 10~* we approach the regime of capillary fingering. In slow capillary fingering the 
local capillary pressures of the menisci becomes equal when all the menisci are stable. This is often referred to as 
capillary equilibrium leading to Peg — Pcf - However, in fast displacement the local capillary pressures are generally 
different. Thus, for a large system with many menisci the average capillary pressure of the front will approach the 
average threshold pressure of the tubes. 

Due to the less viscous defending fluid the pressure gradient at high injection rate is largest at the finger tip closest 
to the upper boundary of the lattice. The meniscus in the uppermost finger tip will therefor more likely continue and 
invade the next tube compared to the menisci lying behind it. Thus, the menisci lying behind will get stuck and the 
fluctuations in the global capillary pressure can only be caused by the moving finger tip. This is seen in Fig. [s] where 
we have plotted Peg (a) together with the capillary pressure of the meniscus which is located in the tube that always 
has the largest flow rate. The data in Fig. ^ are taken from the simulation performed at Ca — 1.1x10"^ given in 
Table |. 

On basis of Figs. |^ and ^ we conclude that the fluctuations in Peg correspond to the local capillary pressure of the 
menisci invading the tubes. Moreover, the average level of Peg is almost equal Pc/. 



5 



B. Stable displacement, M > 1 



In the regime of stable displacement we have run seven simulations spread over different lattices of size between 
25 X 35 and 60 x 60 nodes. See Table ||. The simulations were run with the viscosity ratio M — 1.0x10^, but at 
different capillary numbers. The resulting structure of a simulation performed on a lattice with 60 x 60 nodes is 
visualized in Fig. ^. The pressure across the lattice AP, the global capillary pressure Peg and \/A were calculated 
for all simulations. The results are plotted in Figs. |l^ and |ll| . 

In stable displacements the fluid movements are dominated by the viscous forces in the invading liquid. The viscous 
pressure gradient in the invading phase is found to stabilize the front and a compact pattern with an almost flat front 
between the non- wetting and wetting fluid is generated (see Fig. ^). The stabilized front introduces a length scale in 



the system for large times. This length scale is identified as the saturation width Ws of the front 25 1. 

To save computation time the simulations performed on the lattices of 40 x 40 and 60 x 60 nodes (Fig. |l^) were 
stopped after the width of the front had stabilized. The other simulations (Fig. ^l]), however, were run until the 
invading fluid reached the outlet. 

The pressures in Figs. |l^ and |ll| show a quite different behavior compared to the pressures in Figs. I and 0. The 
pressure across the lattice AP (a) and the global capillary pressure Peg (b) are both found to increase as the more 
viscous fluid is injected into the system. However, the average slope depend very much on the injection rate. At low 
capillary numbers {Ca = 7.5 x 10~^) we reach the regime of capillary fingering and similar to viscous fingering the 
pressures reduce to that describing the capillary fluctuations along the front. 

To explain the rather unexpected behavior at high injection rates we have to discuss the effect of the trapped cluster 
of the defending fluid left behind the front. In stable displacement the driving pressure gradient lies between the inlet 
and the front causing a pressure drop between the top and the bottom of the trapped clusters. At moderate injection 
rates where the clusters stay in place and keep their shapes, the forces due to the pressure drop must be canceled by 
capillary forces acting on the cluster menisci. The sum of those capillary forces contribute to the observed increase 
in Peg. 

The above interpretation of Peg also applies in the regime of viscous fingering. In viscous fingering at high capillary 
numbers there are few clusters. Thus, the extra contribution to Peg from the cluster menisci becomes negligible. When 
the injection rate is reduced clusters develop, but now the pressure gradient across the clusters is small because of the 
low injection rate. Therefore, the contribution to Peg still is negligible. This is in agreement of the observations in the 
previous section where Peg was found to fluctuate around a constant level for both high and low capillary numbers. 

All the displacements listed in Table ||, except that at Ca — 7.5 x lO^"', reached the saturation time tg where the 
front stabilize. In Figs. |l^ and |ll| tg is indicated by a vertical dashed line. If we neglect the fluctuations, AP and Peg 
are approximately linear functions of time for t > t^. 

The linearity in Peg is explained by looking at the generation of the clusters in the system. We notice that for 
large times {t > t^), when the front has saturated with fully developed clusters behind it, the average front position h 
can only depend on the injection rate and the viscosity ratio. Both properties are constant through each simulation, 
resulting in h being linear in time. A plot showing the average front position as a function of time at Ca = 4.6 x 10"'^ 
for the lattice with 60 x 60 nodes is shown in Fig. 

On average, every cluster contributes to Peg with a certain amount causing the global capillary pressure to be 
proportional to the number of clusters behind the front. For large times, the number of clusters increases linearly 
with the average front position causing Peg to be linear in h. Summarized we obtain 

Peg oc h (Xt , t > ts . (12) 

It has to be emphasized that there might be large deviations from the observed linearities when clusters of size 
comparable to the system develops. The argument does not apply when t < t^, either. Then the average change in 
AP and Peg depends on the fractal development of the displacement structure. 

The quantity 1/A, defined in Eq. (||), is calculated by using Eq. (0). In Figs. |l| and |ll] we have plotted Aq/A 
expressed as 

Here Aq is equal to the proportionality factor between Q and AP when only the defending fluid flows trough the 
lattice, i.e. Aq = T,K/^iL. K denotes the absolute permeability of the lattice and L is the length of the system. From 
Figs. |l^ and |l[ we conclude that at high capillary numbers and for t > ts, 1/A is proportional to the injection time. 

The linearity in l/A is related to the displacement structure that develop during injection of non- wetting fluid. Let 
us return to equation Eq. IS) which we can rewrite to obtain 
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= SG(AP-Peg)^ , (14) 



where G = AL/Ti defines the simulated mobiUty of the system. It must be emphasized that this mobihty may be 
hard to measure experimentally and it should not be confused with the conventional effective mobility of the system, 
denoted by Ge// (see below). Eq. ( |l^ looks very similar to Eq. (^) describing the fluid flow through a single tube. If 
we now assume that an equation of the form (^ applies on the whole lattice we interpret G as 

K 

G = . (15) 

Me// 

Here ^e/ / is the effective viscosity of the whole system, ji^f / is not trivial to calculate because of the fractal structure 
of the trapped clusters and front. However, if we denote the normalized average front position as h and assume that 
the viscosity of the region behind the front can be expressed by an effective viscosity /12', Me// becomes 

^leff = (1 - h)^ii + hii2 ■ (16) 

Here 0</i<l. By combining Eqs. (|l6|) and (|l5|), we finally obtain after some reorganization 

{M,-l)h + l] , (17) 



1 11 

G °^ 1 ~ 



where = /i2'/Mi defines the effective viscosity ratio. 

To a first approximation, we have calculated ^2 directly from the displacement structure by taking the sum of the 
fluid viscosities multiplied by their respective macroscopic saturations behind the front. However, the calculated 
was not consistent with (Me — 1) given by the average slope of Aq/A plotted in Figs. |l^ and |ll|. This indicates that 
^2 not only depend on the fluid saturations but also the structure of the liquids behind the front. Thus, microscopic 
properties of the fluid configurations has to be considered in order to find the correct Me- 

We are now in a position to derive a more exact formalism describing the pressure evolution and the changes in the 
effective mobility Geff of the system when clusters develop. Geff corresponds to the ratio of the effective permeability 
and the effective viscosity of the system. If we again look at the global capillary pressure Peg, it can be interpreted 
as a sum of the capillary pressure along the front menisci, Pm/ and the capillary pressure of the cluster menisci, Pmc- 
Furthermore, for large times {t > tg) we can express P,„c — ^mc h giving 

Peg = A^c h + Pmf . (18) 

Here Amc denotes the proportionality factor between Peg and h where < ft- < L is the average front position. In 
the limit of very low injection rate the capillary pressure of the front and cluster menisci will approach capillary 
equilibrium causing A„ic — > and Pm/ — Pcj — Peg- 

We can use the relation for P^g in Eq. (|l^) to deduce a formula for Gg/ / . What we are seeking is an equation for 
the flow rate C/, on the form of Darcy's law 

C/-Ge//(AP-P,„;)l . (19) 

The pressure gradient (AP — Pmf)/L accounts for the capillary pressure due to the front, Pmf- 
We start by inserting Eq. (|l8|) into Eq. (|lj) to get 

C/ = (AP - A„, h - Praf) y , (20) 

where the flow velocity U = Q/T,. Eq. ( pO| ) can be written into the from of Darcy's to obtain 

U ^ (AP - P^f) J , (21) 

Me// + Mmc -tJ 



where 

K 



A,nc h . (22) 



7 



From Eq. (21) we immediately interpret ^mc as the increase in viscosity of the invaded region caused by the cluster 
formation. However, this may just as well be seen as a decrease in the effective permeability behind the front [^ . 
Note that the U dependency in Eq. ( p^ ) only indicates changes in Amc between displacements executed at different 
injection rates. The behavior when the flow rate changes during one displacement is not discussed here. From Eq. 
( pl| ) we finally define the effective mobility of the lattice as 

Geff ^ ^ . (23) 

h^ejf ~r ferric 

There is one important interpretation of Ge// when the average front position has reached the outlet, i.e. h = L. 
Then only the invading fluid flows through the system and G^f / becomes equal to the effective mobility of the invading 
phase. If we insert Eqs. ( [l6| ) and (|2^) into Eq. (^3|) we obtain when h — L 



Geff{h = L) = ^^^ . (24) 
M2 H — 

From this expression follows directly a relation for the relative permeability of the invading phase, fc^i defined as 

kr. = ^-^^ = . (25) 

K 

Note that ^2 = l^eff{h = 1) from Eq. (|l6|) and that ^2 is given by the slope of 1/A in Eq. (pT[). 

The above formalism resulting in Eq. (|23|) , takes into account the capillary pressure of the cluster menisci as well as 
the capillary pressure along the front. However, at moderate injection rates where the clusters stay in place and keep 
their shapes, we have found Ge// directly in the simulations by assigning zero permeability to tubes belonging to the 
trapped clusters. Thus, the clusters will be frozen to their initial positions and in the calculations they are treated as 
additional boundary conditions where fluid cannot flow. Simulations have provided evidence that when the clusters 
are frozen the parameter A in Eq. (H) adjusts such that the simulated G in Eq. (|lj) becomes equal to the calculated 



Geff in Eq. (E3[). Ge// is calculated by using Peg and AP which are plotted in Figs. |lO|and 11, Moreover, with frozen 



clusters the global capillary pressure Peg, reduces to that describing the capillary pressure along the front, Pmf- 



C. Capillary fingering, M = \ 

In the regime of capillary fingering we have run 17 simulations with viscosity matching fluids (M = 1.0) spread 
over six different capillary numbers. The different capillary numbers and the corresponding injection rates are listed 



in Table [IF The lattice size was 40 x 60 nodes for all simulations. Due to long computation time we only did two 
simulations at the lowest capillary number. For all the other capillary numbers, we ran three simulations. 

In capillary fingering the displacement is so slow that the viscous forces are negligible, with the consequence that 
the main force is the capillary one between the two fluids. Only the strength of the threshold pressure in a given 
tube decides whether the invading fluid invades that tube or not. Since the radii of the tubes (which determine the 
threshold pressures) are randomly chosen from a given interval, the non-wetting fluid flows along the path of least 
resistance. 

The displacement structure of one of the simulations at lowest Ga = 4.6 x 10~^ is shown in Fig. IJ. We observe that 
the invading fluid creates a rough front with trapped clusters that appear at all scales between the tube length and 
the maximum width of the front. 

For all simulations AP and Peg were calculated. The result for six of the simulations each at one of the different 
capillary numbers are shown in Fig. |l^. In the figure (a) denotes AP and (b) denotes Peg. Note that Peg has been 
subtracted by lOOOdyn/cm^ to avoid overlap of the curves at low capillary numbers. 

The front was found to stabilize at high injection rates p^ ] and the saturation time tg is indicated by the vertical 
dashed line in Fig. ^ At the lowest capillary number it might be difficult to estimate the saturation time accurately. 
At very low injection rate the width of the front probably approaches an upper cut off equal to the finite size of the 
lattice. 

At high capillary numbers for t > tg, AP and Peg are found to increase linearly as a function of time. This is 
consistent with the result from stable displacement. However, the difference AP — Peg, is constant through each of 
the simulation opposed to what we observed when M > 1. Fig. |l^ shows the difference AP — Peg for one of the 
simulations at Gq = 2.3x10"^. The plot shows the normalized value of AP — Peg subtracted by 1 such that the 
resulting data fall close to zero. The fluctuations appearing in the 9th digit are caused by numerical round off errors 
in the simulations. 
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In the previous section we suggested that the mobiHty G defined in Eq. (|5|), was a function of the effective viscosity 
/ie//. Moreover, we assumed that /ie// depends on the saturation of the invading and the defending fluid and the 
displacement structure. When /ii = /i2 none of these assumptions apply. From Eq. (|l^ the constant pressure 
difference AP — P^g in Fig. |lj implies that the simulated mobility G becomes constant, even with respect to the local 
displacement structure. The effective viscosity reduces to the viscosity of the two liquids, ^^eff = Mi = Ai2 giving 
A = Aq = TjK/iiiL. The difference AP — Peg, however, depend on the injection rate Q. This is observed in Fig. 
where AP — Peg increases for increasing capillary number. 

When /ii = ^2 the pressure difference AP — Peg has a simple interpretation. It corresponds to the viscous pressure 
drop arising when viscous fluids are moving. When the injection rate is reduced the viscous pressure drop becomes 
negligible and we approach the regime of capillary fingering. At low injection rates the pressure gradient across the 
trapped clusters also vanish giving only a small increase in AP and Peg. This is observed in Fig. |l^ at = 4.6 x 10~^ 
where the average increase in Peg becomes quite small and AP ~ Peg. In the same figure, the sudden jumps in 
the pressure function identify the bursts where the invading fluid proceeds abruptly. This corresponds to Haines 
jumps (§J2|. 

The property that the mobility G is constant when the liquids have equal viscosities simplifies the computation of 
the nodal pressures in the lattice. By substituting A with Aq in Eq. we find that the injection rate is given by 



^oAP + B 



(26) 



This equation has only one unknown, the term B, opposed to the original Eq. having two unknowns, both A and 
B. To verify the result when A is replaced by Aq we have compared the solution found from Eq. (26), necessitating 
one solution of the fiow equations, with the one given when Eq. (|^) is solved twice. Not surprisingly, there is excellent 
agreement between these two results. 

The strong evidence that A = Aq is constant in Eq. (|2^) can be deduced from simple considerations of the energy 
dissipation in the system. In analogy with electrical circuits, we define the total energy dissipation W in the system 
as 

W = QAP . (27) 
The total dissipation must equals the sum of the dissipation in every tube a in the lattice. Thus, 

gAP = J2 1'-^P'- ' (28) 



where the summation index a runs over all tubes in the lattice, qa is given by Eq. (^ which we rewrite as {ij a) 

qa = OqAPq + ba ■ (29) 



Here we note that Oq, is proportional to ka/ ^eff, the mobility of the tube a. By inserting Eq. (|5D and ( P9[ ) into Eq. 
(Eq) we get after some reorganization 



AAP + B = ^a„AP('^V+5„^^" 



AP y " AP ■ 

By replacing the local pressure Ap^ in Eq. ( ^0| ) with AP, using Eq. (^ we obtain after some algebra 



(30) 



Q 



H, 



AP 



^ra(2aaHQ + ba 



AP 



(a^Ha + ba) 



(31) 



If we compare the above equation with Eq. (g) we recognize the first summation as A and the second as B. Thus, A 
depends entirely on aa and Fq. As stated earlier, both and F^ are proportional to the mobility of the tubes. The 
mobility of each tube depends on the local fiuid configurations, through the effective viscosity ^eff as defined in Eq. 
(p^). However, when the fluids have equal viscosities we get /^e// = Mi =1^-2- As a consequence aa and Fq, becomes 
constant, which is consistent with the simulation result. 
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IV. DISCUSSIONS 



We have simulated drainage displacements at different injection rates for three different viscosity ratios M — 
1.0 X 10"'^, 1.0 and 1.0 x 10^. The main focus of the work is the study of the temporal evolution of the pressure when 
a non-wetting fluid displaces a wetting fluid in porous media. Moreover, the effect of the trapped clusters on the 
displacement process has been discussed. From the results we clearly see that the capillary forces play an important 
role at both high and low injection rates. 

At high injection rates with M > 1.0 the global capillary pressure, Peg was found to increase as a function of 
the number of trapped clusters behind the front. For large times when the front has saturated. Peg becomes even 
proportion to the average front position h. This lead to a formalism describing the evolution of the effective mobility 
when the non-wetting fluid was injected into the system. Moreover, we showed that this effective mobility could 
be used to estimate the relative permeability of the invading phase when the average front position has reached the 
outlet. 

At moderate injection rates the effective mobility given by Peg, were shown to be equivalent to assigning zero 
permeability to tubes belonging to the trapped clusters. When M ^ 1.0 or at low injection rates the effect of the 
clusters become negligible reducing Peg to describe the local capillary fluctuations of the invading menisci along the 
front. With displacements performed with equal viscosities, AI — 1.0, we found that the difference AP — Peg was 
constant. This was shown to be consistent with the energy dissipation in the system. 

It must be emphasized that the properties we report are only valid for drainage. So far the model is not capable to 
simulate imbibition. The simulations have also been performed on a two-dimensional porous system where clusters 
develop more easily, compared to fluid flow in three-dimensional porous media pst . Moreover, the lattice sizes are 
limited by the computation time and more sophisticated and efficient algorithms have to be developed in order to 
increase the system sizes and thereby improve the above results. Another, and not less important exercise is to 
compare our simulation results with experimental measurements. 
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FIG. 1. Flow in a tube containing a meniscus. 

FIG. 2. Four different fluid arrangements inside one tube. The shaded and the white regions indicate the non-wetting and 
wetting fluid respectively. 

FIG. 3. The motion of the menisci at the nodes, (a) The non-wetting fluid (shaded) reaches the end of the tube (position 

1) and is moved a distance 5 into the neighbor tubes (position 2). (b) The wetting fluid (white) roaches the end of the tubes 
(position 1) and the non-wetting fluid (shaded) retreat to position 2. To conserve the volume of the fluids a proper time is 
recorded due to the small movement 5 in (a) and (b). 

FIG. 4. The displacement structure obtained from a simulation at Ca = 4.6x10"' and M = 1.0. The size of the lattice is 
40 X 60 nodes and the invading non-wetting fluid (black) displaces the defending wetting fluid (gray) from below. Notice the 
rough front between the fluids and the trapped cluster of defending fluid left behind. 

FIG. 5. The displacement structure obtained by a simulation in the regime of viscous fingering on a lattice of 60 x 80 nodes. 
Ca = 4.6 X 10~^ and M = l.Ox 10~^. The invading, non- wetting fluid (black) displaces the defending, wetting fluid (gray) from 
below. 



FIG. 6. AP (a) and Peg (b) plotted as a function of time. Ca = 4.6 x 10"^ and M = l.Ox 10"^. 

FIG. 7. Pcf (a), and Peg (b), at four different capillary numbers for the simulations with M = 1.0 x 10~^ at lattice size 25 x 35 
nodes. The pressures are normalized using the average threshold pressure of the tubes. Note that Peg has been subtracted by 
1000 dyn/cm^ before it was normalized to avoid overlap between the two curves at low capillary numbers. 

FIG. 8. Peg (a) and the capillary pressure of the meniscus traveling with the highest velocity (b) at Co = 1.1x10"^ and 
M = l.Ox 10~^ in the time interval between 7.0 and 12.0s. 



FIG. 9. The displacement structure obtained of a simulation in the regime of stable displacement on a lattice of 60 x 60 
nodes. Ca = 4.6 x 10~^ and M = l.Ox 10^. The invading, non-wetting fluid (black) displaces the defending, wetting fluid (gray) 
from below. 

FIG. 10. AP (a), Peg (b) and Aq/A (c) as functions of time for the simulations with M = l.Ox 10^ at lattice size 60 x 60 
(to the left) and 40 x 40 nodes (to the right). The vertical dashed lines indicate the saturation time where the front stabilize. 

FIG. 11. AP (a), Peg (b) and Ao/A (c) as functions of time for the five simulations with M = l.Ox 10^ at lattice size 25 x 35 
nodes. Note that for Ca = 4.2x10"^, 2.2x10"^ and 7.5x10""' Peg has been subtracted by lOOOdyn/cm^ to avoid overlap of 
the curves. The vertical dashed lines indicate the saturation time where the front stabilize. 
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FIG. 12. The average front position as a function of time at Ca ~ 4.6 x 10 ^ and M — 1.0 x 10^ for the lattice of 60 x 60 
nodes. 



FIG. 13. AP (a) and Peg (b) for six simulations with AI = 1.0 each at one of the capillary number listed in Table [II 
The dashed line indicates the saturation time ts when the front has reached the saturation width w, 
subtracted by 1000 dyn/cm^ to avoid the curves from overlapping at low capillary numbers. 



FIG. 14. The normalized difference AP - Peg subtracted by 1 at Ca = 2.3x10"^ and M = 1.0. 



TABLE I. The lattice size and the values for the injection rate and the capillary number when M = l.Ox 10 



Size 


Injection rate 


Ca 


(nodes) 


(cm'^/min) 




60 x 80 


1.5 


4.6x10"^ 


25 X 35 


1.4 


1.1x10"^ 


25 X 35 


0.98 


7.1x10"^ 


25 X 35 


0.62 


4.7x10"^ 


25 X 35 


0.50 


3.6x10"^ 


25 X 35 


0.099 


7.2x10"* 


25 X 35 


0.049 


3.5x10"* 



TABLE II. The lattice size and the values for the injection rate and the capillary number when M — 1.0x10^. 



Size 


Injection rate 


Ca 


(nodes) 


(cm^/min) 




60 X 60 


1.5 


4.6x10"'' 


40 X 40 


1.0 


4.6x10"^ 


25 X 35 


2.5 


1.8x10"^ 


25 X 35 


1.3 


9.5x10"^ 


25 X 35 


0.57 


4.2x10"^ 


25 X 35 


0.29 


2.2x10"^ 


25 X 35 


0.10 


7.5x10"* 



TABLE III. The values for the injection rate and the capillary number when M = 1.0. The lattice is 40 x 60 nodes. 



Runs 


Injection rate 
(cm^/min) 


Ca 


3 


10 


2.3x10"^ 


3 


4.0 


9.2x10"* 


3 


2.0 


4.6x10"* 


3 


1.0 


2.3x10"* 


3 


0.40 


9.2x10"^ 


2 


0.20 


4.6x10"^ 



12 



(b) (c) (d) 



(a) 



(b) 



0.50 




50 100 150 

Time (s) 




0.0 2.5 5.0 7.5 10.0 12.5 
Time (s) 




0.0 20.0 40.0 60.0 80.0 
Time (s) 




100 200 300 
Time (s) 



400 



c 



3 



73 



73 



3 



5.0 
4.0 
3.0 
2.0 
1.0 
0.0 
-1.0 



■C =9.2* 10"* ! 


(a) . 
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